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Abstract 

We revisit a classical perturbative approach to the Hamiltonian related to the motions 
of Trojan bodies, in the framework of the Planar Circular Restricted Three-Body Problem 
(PCRTBP), by introducing a number of key new ideas in the formulation. In some sense, 
we adapt the approach of m to the context of the normal form theory and its modern 
techniques. First, we make use of Delaunay variables for a physically accurate representation 
of the system. Therefore, we introduce a novel manipulation of the variables so as to respect 
the natural behavior of the model. We develop a normalization procedure over the fast 
angle which exploits the fact that singularities in this model are essentially related to the 
slow angle. Thus, we produce a new normal form, i.e. an integrable approximation to the 
Hamiltonian. We emphasize some practical examples of the applicability of our normalizing 
scheme, e.g. the estimation of the stable libration region. Finally, we compare the level 
curves produced by our normal form with surfaces of section provided by the integration 
of the non-normalized Hamiltonian, with very good agreement. Further precision tests are 
also provided. In addition, we give a step-by-step description of the algorithm, allowing for 
extensions to more complicated models. 


1 Introduction 

Series expansions in terms of small physical parameters are a common way of dealing with 
dynamical models in Celestial Mechanics. One case where this approach has been extensively 
used is the problem of Trojan motion. The so-called Trojan stability problem, i.e. the study of 
the long-term dynamics for massless particles in the neighborhood of the equilateral Lagrangian 
points, has been studied both analytically and numerically in several context^. 

* Key words and phrases: Restricted three-body problem, normal forms, Hamiltonian perturbation theory, 
averaging. Celestial Mechanics 

^For an introduction to Trojan dynamics, see, e.g., [9] 
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From the numerical point of view, quite complete models (including extra planets perturb¬ 
ations and/or 3-dimensional configurations) are possible to consider. Such experiments are in 
general based on accurate long integrations, done through high precision integrators. A clear 
outcome is the determination of size and shape of the stable domain around the equilibrium 
points. Secondary resonances embedded within the domain of the 1:1 Mean Motion Resonance 
(MMR) are of particular importance in the stability issue. For example, in [22] and |23|, the 
importance of the resonance web is highlighted for what concerns the process of de-population 
of the Trojan domain, in the case of Jupiter’s asteroids. In [6], the authors study the orbital 
stability of the only discovered Trojan asteroid of the Earth, i.e. 2010 TK7, and determine the 
stability region around L4 and L5 in the Sun-Earth system. 

Besides pure numerical investigations, several works have also explored the applicability 
of analytical methods in the problem of Trojan stability. In particular, the Trojan problem 
has served as a classical model for testing methods based on the construction of a so-called 
Hamiltonian normal form approach. 

Semi-analytic methods share a common structure: each method is based on an explicit 
algorithm, possible, as a rule, to translate into a programming code, which computes the ex¬ 
pansion of a suitable normal form that provides some special solutions. Therefore, a normal 
form is a good local approximation to the complete Hamiltonian, used to reproduce the orbits 
of particular objects. 

Normal forms have been used in several cases for computing the size and shape of the stability 
domain. From the point of view of Nekhoroshev stability estimates [I1],[T9|, examples of this 
computation are in ID, ii3i, m, m, m- On the other hand, in |10|, the stability estimation is 
induced by applying the Kolmogorov-Arnold-Moser (KAM) theory. But all these approaches, 
while introducing novel and promising features, tend to fail when it comes to reproduce the 
numerical results about the size and shape of the whole stability region. In general, they are 
able to justify stability for just a small subset of Trojan orbits, located relatively close to the 
equilibrium point. 

In the present work, we develop a new normal form approach to the classical PCRTBP by 
combining three main ideas. First, we make an accurate choice of variables for the representation 
of the system. Second, we explicitly take benefit of the existence of a slow and a fast degree of 
freedom, dealing with them differently in our normal form construction. Finally, we make use of 
Lie series techniques for the averaging over the fast angle; in fact, as shown below, a key element 
of our method is that such a treatment makes possible to overcome the only true singularity 
that the Trojan motion contains, namely the eventual close encounters with the primary. 

We motivate these ideas as follows. 

Many of previous attemps describe the system in terms of cartesian coordinates, which 
however do not capture properly the physieal configuration of the model. Considering canonical 
variables well adapted to the system helps to improve the accuracy of the normal form and the 
estimation of stability domains. To this end, in the present paper, we adopt (and show the 
usefulness of) a modihed version of Delaunay-like coordinates introduced long ago in [llj . 

Regarding the other two ideas, an important point is that the only real physical singularity 
in the 1:1 MMR region is due to possible collisions / close encounters of the Trojan body with 
the primary. In our setting of modified Delaunay-like variables (see definition Hj), this singularity 
takes place at A = 0, where A is the synodic mean longitude of the massless body. Thus, any 
polynomial expansion around the equilibrium point, exhibits a bad convergence behavior for 
orbits approaching this singularity. Nevertheless, in this work we show this problem can be 
overcome in a rather simple way. This is based on the following remarks. 
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On one hand, as noted already, an inspection of the main terms of the Hamiltonian for a 
Trojan body, indicates that the motion is ruled by a fast and a slow frenquency. In Delaunay 
variables, these two dynamics are represented by two independent pairs of canonical coordinates. 
However, the long term behavior of the orbits depends essentially only on the slow degree of 
freedom which, in our variables, corresponds to the canonical pair including the angle A. In other 
words, an appropriate normal form construction that aims to study the long term dynamics 
involves averaging over the fast angle only. This produces an integrable Hamiltonian of two 
degrees of freedom in which the fast angle is ignorable. 

On the other hand, such an averaging implies that one has to solve a so-called homological 
equation in which A plays no role. As shown in Section 2 below, due to this property, one can 
retain in the Hamiltonian a complicated functional dependence on A, other than trigonometric 
or simply polynomial. In our case this dependence turns to be of the form of powers of the 
quantity /3(A) = , ^ Thus, using this technique allows to greatly extend the convergence 

domain of the final normal form. 

Finally, we develop a new normalization algorithm for the computation of the integrable 
approximation, by using the Lie series formalism, adapting to a modern way the technique 
described in m- With this integrable normal form, we can approximate Trojan orbits having 
either tadpole or horseshoe shapes, even in cases where distances from the equilateral point 
become large and where previous approaches tend to fail. 

There are several different examples where such a normalizing scheme can be applied. The 
most evident corresponds to the use of the scheme for the estimation of the stability domain. 
Such a computation, in the same direction as the references mentioned before, aim to estim¬ 
ate the effects induced by the remainder TZ (in later Eq. ITT]) of the normal form produced by 
the algorithm. According to our results, our novel normalization may radically improve the 
results regarding the size of the domain of stability around the equilateral Lagrangian points, 
understimated so far. Other examples inherent directly to the normal form can be mentioned. 
For instance, in EH, we used the normalized coordinates for the design and optimization of 
maneuvers aiming to transfer a spacecraft into the tadpole region. Moreover, the numerical ex¬ 
periments described in the present work highlight that our method allows to clearly differentiate 
the tadpole region from the horseshoe region. In some cases, when the mass ratio between the 
two primary bodies is very small and consequently also the chaotic regions are small, such a 
computation gives a first order estimation of the stable tadpole domain. On the other hand, this 
kind of perturbative approach is motivating for problems of diverse nature. For instance, in [3], 
they make use of the Relegation algorithm, that shares a similar structure with our normalizing 
scheme, to compute some particular orbits around an irregularly shaped asteroid. 

This paper is structured as follows: in Section[21 we describe the expansion and normalization 
scheme, applied to the case of the PCRTBP Hamiltonian, in the 1:1 MMR region; in Section [3l 
we provide different tests for a suitable accuracy verification of the normal form, involving 
comparisons with the numerical integration of the full problem; in Section 01 we summarize 
the work and outline future applications of this method. Furthermore, Appendix formally 
presents the algorithm in such a way that it can be adapted also to different models, for example, 
the Elliptic Restricted Three Body Problem (ERTBP) or the Restricted MultiPlanet Problem 
(RMPP), described in [20]. 
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2 Construction of the integrable approximation 


2.1 Initial settings 

In heliocentric canonical variables (p, r), the Hamiltonian of the PCRTBP can be written as: 



GM 

r 




( 1 ) 


where M and m' are the masses of the larger and smaller primary, respectively, r is the he¬ 
liocentric position vector of the test particle, r' the one corresponding to the second primary, 
p = IIpII, r = ||r||, r' = ||r'|| and A = ||r — r'||. 

In what follows, a, e, M and w symbolize the major semi-axis, the eccentricity, the mean 
anomaly and the longitude of the perihelion (primed quantities correspond to the primary). We 
set the unit of length as a' = 1, and the unit of time such that the mean motion of the primary 
is equal to n = 1. This implies G{M = 1. The unit of mass is set so as to G = 1. Defining 
the mass parameter as p = m!, the Hamiltonian in the above units takes the form: 

1 

H = ^ - ^^F , ( 2 ) 

1 r 


where the so-called disturbing function jjlF is such that 

Y A r'3 r) 


(3) 


Including the small keplerian correction /x/r in the disturbing function allows to define action- 
angle variables with values independent of the mass parameter p. Inspired by [20] and references 
therein, we introduce modified Delaunay action-angle variables 


G = \/ a(l — e^) , X = M + zu — M' — w' , 

r = ^(i- y/i - e2^ , £ = M . 

According to the definitions of the orbital elements above, A is the synodic mean longitude. In 
order to remove the fictitious singularity of the Delaunay-like variables when e = 0 (i.e., for 
circular orbits), it is convenient to introduce canonical coordinates (p, A, p), similar to Poincare 
coordinates for the planar Keplerian problem. Thus, let us define 

P = G — I , A = A, 

^ = \/^ cos i , p = sin £ , 


where the values of p, ^ and p are significantly small in a region surrounding the triangular 
Lagrangian points, for instance, in the case of tadpole or horseshoe orbits (since a ~ 1 and 
e > 0). Let us recall that those variables have been introduced in [11] to study the PCRTBP. 

Before starting the construction of the normal form, it is necessary to express Hamilto¬ 
nian ([2]), in particular the disturbing function F in ([3]), in terms of Poincare-Delaunw-like 
coordinates {p,^,X,p). Here, we limit ourselves to sketch this preliminary procedurqj. Es¬ 
sentially, as a first step, the terms 1/r, and r • r' (appearing in the Hamiltonian) must be 

^An exhaustive description of those expansions is out of the scopes of this publication. All the details are 
deferred to Paez, R.I., Ph.D. Thesis (2015), that is in preparation and is available on request to the author. 
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expanded with respect to orbital elements a, e, M and A, following e.g. §6 in HZ]. Afterwards, 
the eccentricity e is replaced by its expansion in power series of the parameter s/TT/G. By 
using the equations a = {G + B)^, G = 1 + p, F = (^^ + p^)/2, we can then express 1/r, and 
r • r' as functions of the canonical variables ([5|) . 

As explained in the Introduction, the particular normalizing scheme we use allows to keep 
a non polynomial nor trigonometric functional dependence with A. By making use of this idea, 
we keep powers of /3(A) = l/y/2 — 2 cos A, in the expansions of the term ^ of F, which describes 
the main part of the inverse of the distance from the primary (for values of the major semi-axis 
a ~ 1 and small eccentricity). Thus, the Hamiltonian ruling the motion of the third body takes 
the following formH 

'H{p,C,\v)= “2 —) - 1-P +/^(l + cos(A)-/3(A)) 

^ j=o ^ ^ 

oo (6) 

l=l mi+m2 k]_+k2<l 
+m3=Z i<2/+l 

where bmi,m 2 ,m 3 ,ki,k 2 ,j rational numbers. Actually, we just produce a truncated expansion of 
formula ([6]) up to a finite polynomial degree in p, ^ and p, by using Mathematica. By inspection, 
we check that both properties ki + k 2 < I and j < 21 + 1 are satisfied in our initial expansion, 
and we just conjecture that they hold true also to any following polynomial degree in p, ^ and 
p (being such a proof beyond of the scopes of this work). 


2.2 Algorithm constructing the normal form averaged over the fast angle 

Let us focus on the first main terms of the Keplerian part: 


3 e+d^ 

2^ 2 


P + 


+ p'^ 


+ 



5(p + r)^ + ... , 


(7) 


where we refer to the dynamics of the new canonical coordinates (^, p) by retaining the old 
action-angle pair (F,!"), that allow us to sketch our strategy in a simpler way. Indeed, the 
previous formula shows that the angular velocities have different order of magnitude 


A 



t = 


dn 


~ I, 


( 8 ) 


because /r <C 1 and the values of the actions p and F are small in a region surrounding the 
Lagrangian points. Thus, A can be seen as a slow angle and ^ as a fast angle. This motivates 
to average the Hamiltonian over the fast angle (see, e.g.,§ 52 of 121), in order to focus mainly on 
the secular evolution of the system. Therefore, we want remove (most of) the terms depending 
on the fast angle £, that can be conveniently done in the setting of the variables {p,^,\,p), by 
performing a sequence of canonical transformations. In the following, this strategy is translated 
in an explicit algorithm using Lie series, by adapting to the present context an approach that 
has been fruitfully applied to various problems in Celestial Mechanics (for instance, for locating 
the elliptic lower-dimensional tori in a rather realistic four-body planetary model in [21], or 


®The first term in (El) comes from the expansion of the Keplerian part, i.e. 
independent of p, in terms of C V and p, that corresponds to /C = — 


2(p+i+i^) 


the part of the Hamiltonian 

-p- 1. 
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studying the secular behavior of the orbital elements of extrasolar planets in m)- By applying 
such a procedure to our Hamiltonian model, we construct an average normal form, satisfying 
two important properties: (I) it provides an accurate approximation of the starting model; (II) it 
only depends on the actions and one of the angles and, therefore, it is integrable. 

Our final normal form is produced by an algorithm dealing with a sequence of Hamiltonians, 
whose the expansion can be conveniently described after introducing the following 

Definition 2.1 A generic function g = g{p, A, r]) belongs to the class Vi^s , if its expansion is 
of the type: 

E E (cosA)^HsmA)^^ (/3(A))^' , 

2rai+rrL2+m2,=l fci+fc2<^+4s—3 
j<2/+7s-6 


whcTC T6(ll COcfflciCTltS. 

Let ri and r 2 be two integer counters, running in the intervals [1, i?i] and [0, R 2 ], respectively, 
being Ri , R 2 E N fixed numbers. At each (ri,r 2 )-th step, our algorithm introduces a new 
Hamiltonian guch that 


Z>4 




ri —1 / i?2 


2 

2 I ^2 


(s)/ ^ 


1>R2 


s r{ri,r2;s) 


+ X] X] ' [ P 

s=l \l =0 

^2 y-2 2 

An)/ g +h ^ ST . ,ri f{ri,r2-,ri) 


+ EE‘A’ 


1=0 


l>r2 


(/3,C,A,r/) 


(9) 




s>ri l>0 


where E Pip V I > 4, e Pi^s V 0 < I < R 2 , 1 < s < ri, E Pi^n V 0 < I < r 2 , 
jhi,)- 2 ,n) ^ \/ I > r 2 , e \/ I > R 2 , I < s < ri and V / > 0, s > ri. From Q, 

we emphasize that 


• The splitting of the Hamiltonian in sub-functions belonging to different sets Vi^s basically 
gather all the terms with the same order of magnitude and total degree 1/2 (that can 
be semi-odd) in the actions p and F. This is made in order to develop a normalization 
procedure, exploiting the existence of natural small parameters: p and the values of the 
pair of actions (p, F). 

• All the terms Z^^'^ and appearing in equation ([9]) are made by expansions including 

a finite number of monomials of the type described in Definition 12.11 

• At the beginning of our algorithm, we can set = R, because the expansion ([B]) of 

the initial Hamiltonian R can be expressed as in equation ([9|) . 

Our algorithm requires just R 1 R 2 normalization steps, which are performed by constructing 
the finite sequence of Hamiltonians 

= R, ... , ^ h{R 2 ,Ri) . 
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They are defined so that V 1 < ri < Ri and the (ri,r 2 )-th normaliz¬ 

ation step is performed by a canonical transformation. This recursively introduces the new 
Hamiltonian in such a way that 

^ , ( 10 ) 

with a generating function = Xr^ 2 ^\p^ d) that is determined by solving a so-called “ho¬ 
mological equation”, while exp (£^) • = Ylj>o ji^x ' denotes nothing but the Lie series operator. 
The Lie derivative C^g = {g,x} is such that {•,•} is the classical Poisson bracket, with g a 
generic function defined on the phase space and x any generating function (see, e.g., [12] for an 
introduction to canonical transformations expressed by Lie series in the context of the Hamilto¬ 
nian perturbation theory). All the recursive formulas, which determine the terms of type Z and 
/ appearing in Q, are reported in Appendix lAl Let us stress that, after each transformation, 
in the present subsection we do not change the name of the canonical variables in order to 
simplify the notation. We emphasize here that the new generating function introduced at the 
generic (ri,r 2 )-th step, namely is determined so as to remove from the main 

perturbing terrr@ its subpart that is not in normal form. 

We rewrite the final Hamiltonian in such a way to distinguish the normal form part from 
the rest, as follows: 


A, g) = (p, (^2 + 7?2)/2, A) + C, A, g) , (11) 

where all the averaged terms of type Z are gathered into the integrable part while the 

others contribute to the remainder ^,(^^’^ 2 ) ^ Here, our algorithm is described at a purely formal 
level in the sense that the problem of the analytic convergence of the series on some domains is 
not considered. However, it is natural to expect that our procedure defines diverging series into 
the limit of 00 , because a non-integrable Hamiltonian cannot be transformed in an 

integrable one, on any open domain. In principle, the values of both integer parameters Ri and 
i ?2 should be carefully chosen in such a way to reduce the size of Rjy^^T 2 ) g^g ixmch gg possible. 
In practice, we simply hxed the values of Ri and R 2 according to the computational resources, 
in order to deal with the application described in the following section [3j 

2.3 Numerical computation of the flow induced by the integrable approxim¬ 
ation 

Lie series induce canonical transformations in a Hamiltonian framework; this fundamental fea¬ 
ture allows us to design a numerical integration method, by using both the normal form previ¬ 
ously discussed and the corresponding canonical coordinates. Let us denote with (p(^i>'’ 2 )^ 
the set of canonical coordinates related to the (ri,r 2 )-th normalization step. By applying the 
so-called ’exchange theorem’ (see m), we have that 

fp{ri,r2) ^p(ri,r2) ^ \Pl,r2) ^ .^ 2 )^ — ff(ri,r2-l) ^p{ri,r2) ^p(ri,r2) ^ ^(ri,r2) ^ ^(-^1 ,^ 2 ) ^ p(ri ,r2)^'^ ^ 

(12) 

^Let us recall that the size of fT G Tr 2,3 is expected to decrease when the indexes s or r 2 are 

increased, because the values of P a-ud are assumed to be small 
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where the variables related to the previous step, i.e. i)^^(n .’'2 i))^ are 

given as 


(^(n,r2)('^(ri,r2)^^(ri,r2)^^(ri,r2)^^(ri,r2)^ = g^p f/i ^ Vphl.^2), .^2), .’'2)\ _ (13^ 

\ M Xro / 


According to this, Lie series must be applied separatedly to each variable, in order to properly 
define all the coordinates for the canonical transformation Thus, the whole normaliza¬ 

tion procedure is described by the canonical transformation 

^ (^(1,1) o _ o (^(h« 2 ) o _ _ _ o ipiRlR . . . o _ (-14) 


Such a composition of all the intermediate changes of variables can be used for providing the 
following semi-analytical scheme to integrate the equations of motion: 


^(^(-Ri,-R2)^ 




2(«1.«2) 


p(0’0)(t),^(0’0)(t),A(0’0)(t),ry(0’0)(t)) 


^(fll,i?2) 


^(iJl ,iJ2) ^ ^(Kl ,R2) ^ ,K2) ^ ^(Rl ,iJ2) 

(15) 


where <1>^ denotes the flow induced on the canonical coordinates by the generic Hamiltonian 
fC for an interval of time equal to t. Let us emphasize that the above integration scheme 
provides just an approximate solution. From an ideal point of view (i.e., if all the expansions 
were performed without errors and truncations), formula (|15l) would be exact if the normal 
form part 2 ^(^i’^ 2 ) would correspond to the complete Hamiltonian moreover, Z{RiR2) 

is integrable and its flow is easy to comput^, reasons why using 2 ^(^i’'^ 2 ) becomes valuable. 
According to equation m, the approximate solution provided by the scheme (1151) is as more 
accurate as smaller the perturbing part 7 ^(^i’^ 2 ) jg with respect to 2 ^(^i>^ 2 )_ 

A key point concerns the structure of our expansions. If we truncate the r.h.s. of equation ([9]) 
up to a finite order of magnitude in pL and a fixed total degree //2 in the actions p and F = 
(^^-|-?7^)/2, since each term of type Z and / is included in a corresponding class of functions Vps , 
the truncated expansion of include a finite number of monomials. The same applies 

also to the truncated expansions of both the final Hamiltonian ^nd the canonical 

transformation Therefore, after introducing the truncation rules, the normalization 

algorithm require a finite total number of operations, to compute the elements necessary to 
implement the whole integration scheme m- Thus, it can be translated in a programming 
code. 


®In order to explicitly describe the solutions of the equation of motions for the normal form 
it is convenient to introduce the temporary action-angle variables such that _ 

V2r(^i’^2) cosand where is a constant of motion for the 

normal form _ 2,^Ri,R2) (^p{Ri,R2) ^ gy considering as a fixed parameter and 

using the standard quadrature method for conservative systems with 1 d.o.f., one can compute and 

\(Ri,R2 )at any time t. The same can be done for the evolution of by evaluating the integral corres¬ 
ponding to the differential equation practical purposes, the application of the classical 

quadrature method can be replaced by any numerical integrator that is precise enough. Finally, the values of 
f(Ri,R2)u\ and 77)^1’^2) b) can be directly calculated from those of the corresponding action-angle variables, that 
are and 
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In order to concretely apply our semi-analytical scheme, by using Mathematica, we compute 
the truncations up to the terms 0{p^) and to the fifth total degree with respect to the square 
roots of the actions p and -|- ?/^)/2 (i.e. Ri = 3, R 2 = 5, in order to obtain ~ 

Ci3,5) 

inverse). At the end of its execution, that code writes, on an external file, C 
functions which are able to compute the changes of coordinates induced by and by its 

inverse. Moreover, it provides also another external file, where the integrable equations of motion 
related to the Hamiltonian are written according to the syntax of the Tayloi| software 

package. A basic use of the Linux shell scripting allowed us to automatically complement all 
these parts of the computational procedure, so as to produce a new C program as an output. 
This final code is able to numerically integrate the equations of motion via the scheme (I15p , where 
the flow 5 ) is computed by using the Taylor method (based on the automatic differentiation 
technique, see m and [1]). The truncation rules are fixed in such a way to execute all the 
computations in a reasonable amount of CPU-time. 

3 Accuracy control for the normal form 

An averaged model approximating a certain problem is powerful to the extent that it reproduces 
the main features of the original system, (in this case the Planar Circular Restricted Three Body 
problem). It is widely proved that PCRTBP is barely a very simplistic representation for the 
Trojan motion, but its basic dynamics allows us to test in a simple way the normalization 
method. Furthermore, according to the nature of the method, much more complex (i.e., non- 
planar, eccentric, including additional planets) models can be treated without any substantial 
change in the averaging scheme. 

When seen from the most commonly used reference system (i.e. a synodic rotating frame 
with origin in the barycenter), motions around the equilateral Lagrangian points are represented 
by the so-called tadpole or horseshoe orbits (see §3.9 of m)- In our set of variables, these 
orbits are characterized by large variations of the angle A. So, as starting, a suitable averaged 
approximation has to be able to reproduce these variations correctly. In particular, this should 
apply even for the challenging case of bodies whose orbits lay very close to the border of the 
stability domain. As second goal, the averaged Hamiltonian should be able also to distinguish 
between tadpole orbits (around just one equilateral equilibrium point) and horseshoe orbits 
(around both points). In Section [2l we have provided an integrable averaged Hamiltonian 
Z{RiT2) 

(see El), that explicitly describes the behavior of the degree of freedom related to the 
canonical pair of variables (p,A). In order to compare the original problem with our averaged 
version, we develop two different tests. 

3.1 Numerical surfaces of section vs. semi-analytical level curves 

The first test consists of a graphical comparison between the orbits provided by the complete 
Hamiltonian and those provided by normal form both cases, the initial conditions 

are fictitious, but each set is derived from the catalogued position of a real generating Trojan 
body. Since we are just interested in the behavior of the slow degree of freedom, for the complete 
problem (originally 4D), we just retain the evolution of the variables A and p by means of 
isoenergetic surfaces of section, which gives a 2D representation. In the two cases presented 

®Taylor is an automatic translator producing a C function acting as a time-stepper specific for a given 
ordinary differential equation (by means of the Taylor method). It is publicly available at the following website: 
http: //www.maia.ub. es/~angel/soft .html 
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below, the generating bodies are 2010 TK 7 , Trojan asteroid of the Sun-Earth system and the 
asteroid 1872 Helenos of the Sun-Jupiter system. 

From a catalogue, we obtain the coordinates (rotated to the plane of primaries) of each gener¬ 
ating body for a certain epoch and we convert them to our canonical coordinates {pgb, \gb, Cgb, flgb)- 
This initial condition provides the Jacobi constant CJ^^^ for the body, after which it is possible 
to produce an isoenergetic surface of section. We generate a set of 10 new orbits by keeping 
fixed the initial values for the variables p = pgb and rj = pgb, scanning a range of values for 
A and obtaining ^ in such a way that Cj{p, X,^,r]) = (isoenergetic orbits). These initial 
conditions are numerically integrated for a short time, up to the variables accomplish the rela¬ 
tion M{^,r]) = 0, where M corresponds to the mean anomaly. We call this new set Sgb, and it 
generates both the surface of section and the level curves. 

3.1.1 Computation of surfaces of section 

Starting from Newtonian formulation for the PCRTBP, we derive the equations of motion for 
the third body in cartesian coordinates in the heliocentric systerr 0 

X — Xp ^ Xp 

_ y-yp _^ yp 

where x, y, Vx, Vy correspond to cartesian positions and velocities of the massless body, and xp, 
yp give the instantaneous position of the planet. We translate initial conditions of the set Sgb to 
cartesian coordinates and we integrate them with a Runge-Kutta 7 - order integrator, along 
1500 periods of primaries, with time-step equal to 27r/100. During this integration we collect 
the points contained in the pericentric surface of section, which in our case is represented by the 
condition 

= 0 or equivalently p = 0 , (17) 

and provide about 1500 points per orbit. The output data is again translated to Delaunay 
variables and suitable to be compared with the results of the averaged Hamiltonian. 





y =Vy , 


(1 - y)y 




3 . 1.2 Computatiou of level curves 

Through Hamilton’s equations, we can easily derive the equations of motion for the canonical 
variables p and A, 

Q2,^Ri,R2) . Q2,^Ri,R2) 

^ dX ^ dp 

Since is given by a series expansion of the variables p, T, cos A, sin A, /3(A) and the small 

parameter p, the equations of motion inherit the same structure. Through a suitable storing 
of the coefficients of these series and management of the equations of motion implemented in 
Mathematica (as explained in last paragraph of Subsection 12.31) . we compute the evolution of 
the orbits according to the averaged Hamiltonian. Every initial condition of the set Sgb is first 

^Although cartesian coordinates are not canonical conjugated variables, they provided the simplest closed form 
for these equations of motion, and for this reason they are still widely used for numerical experiments. 
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Figure 1: Comparison between the level curves produced by the averaged Hamiltonian in red 
(light gray) and the points of the surface of section for the complete problem in blue (dark gray), 
for the Sun-Earth problem (left panel) and Sun-Jupiter problem (right panel). In the Sun-Earth 
case, the generating body is the Earth Trojan 2010 TKy. In the Sun-Jupiter case, the generating 
body is Trojan asteroid 1872 Helenos. See text for more details. 


converted to normalized variables and then integrated up to collecting about 2000 points, keeping 
the relative energy error of the integration smaller than 10“^^. Such a numerical integration is an 
efficient way to compute the level curves for the integrable normal form corresponding 

to the values F = (^^ + rj^)l2 and ,\) (in normalized variables). 

We complement every point of a level curve with the values ^ = \/^ and rj = 0 (equivalent 
to M = 0). Let us note here that the condition M = 0 in the normalized coordinates does not 
correspond exactly to the surface of section M = 0 in the original variables. However, since 
the change of coordinates (that gives the values of the non~normalized variables in the 

scheme (1151) 1 is, by construction, a near-to-identity canonical transformation, we assume that 
the conditions for each surface of section do not differ too much. Finally, via we back- 

transform all the points of a level curve in the original variables and we graphically compare 
them with the corresponding numerical surface of section. 

3.1.3 Examples and results 

We choose two systems with very different values of the mass parameter for a better contrast 
in the results of the test. The first case is provided by the PCRTBP approach to the Sun- 
Earth system, which is defined by a mass parameter p = 0.30003 x 10“®. The generating body 
chosen for this system is the Earth only Trojan 2010 TK 7 , for which we obtain its coordinates 
from Jet Propulsion Laboratories JPL Ephemerides Servicqj, at epoch 2456987.5 JD (2014-Nov- 
26). We translate them to our set of canonical variables, resulting ptk 7 = —1.8401447 x 10“^, 
^TK 7 = 3.5736334, riTK 7 = 0.1152511 and ^tk 7 = —0.1530054. For the second case, we 
choose the Sun-Jupiter system, defined by the mass parameter p = 0.953855 x 10“^. The 
generating body in this case is the Trojan asteroid 1872 Helenos, which belongs to the Trojan 
camp around L5 in such a system. The set of initial conditions for Helenos were obtained 
from Howell Catalogu^, at 2452600.5 JD (2002-Oct-22), and after being translated to canonical 

^http://ssd.jpl.nasa.gov/?ephemerides 

® http: //www.naic . edu/~nolaii/astorb .html 
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Delaunay variables, they read/9i872 = —0.3836735x10“^, A 1872 = 5.6716748, r /1872 = —0.0154266 
and ^1872 = —0.1104177. 

Figured] shows the comparison between the surface of section and the level curves computed 
for Sun-Earth system (left panel), and Sun-Jupiter system (right panel). For both cases, the 
points of the surface of section are represented in blue (dark gray) and the curves produced with 
the averaged Hamiltonian are in red (light gray). For the case of Sun-Earth system the agreement 
between the two representations is excellent. According to the milestones we dehne at the 
beginning of this section, the averaged Hamiltonian reproduces accurately the large variations 
of A. In particular, it is perfectly able to distinguish between orbits belonging to the tadpole or 
to the horseshoe region. In the case Sun-Jupiter system, for which the mass parameter value is 
3 orders of magnitude larger, there is a substantial presence of chaos. Nevertheless, the averaged 
Hamiltonian is able to locate any tadpole orbit provided by the complete Hamiltonian, even in 
cases when the motion is trapped into secondary resonances, and its validity is also good for 
orbits close to the border of the stable region. 


3.2 Computation of quasi-actions 


So far in the literature, one the most successful attemps to construct a normal form for testing 
the stability of Trojans in the context of the Hamiltonian formalism is discussed in [TO]. However, 
of the 34 initial conditions they used for the system Sun-Jupiter, 4 presented orbits that were 
highly chaotic after being translated according their rotations into the planar CRTBP, while the 
Kolmogorov normalization algorithm defined in that work did not work properly for other 7 cases. 
Here, we revisit the fictitious initial conditions they considered for these latter seven asteroids 
(1868 Thersites, 1872 Helenos, 2146 Stentor, 2207 Antenor, 2363 Cebriones, 2674 Pandarus and 
2759 Idomeneus). Since they conform a set of coordinates that either lay very close to the border 
of stability or show an anomalous behavior (with respect to the expected tadpole orbit), they 
provide a natural harder test in a more quantitative way. 

In the previous subsection, we show that the evolution of the orbits given by emu¬ 

lates correctly in many cases the evolution under the original non-normalized Hamiltonian, by 
means of graphical comparisons between level curves and surfaces of section. This normal form 
2 ^(Ri,R 2 ) contains two different actions or integrals of motion. One, obtained by construction 
and through the normalization, is given by P. The other one, not explicitly obtained, is due to 
the fact that, after the reduction of P, the normal form bear just 1 d.o.f., i.e. it is integrable. 
This second constant of motion is associated with the area enclosed by the level curve computed 
under the integration of and therefore, provide another quantity to be checked. In 

order to do so, the computation of the orbits is done as explained in subsections I3.1.11[TTT31 
but for just one initial condition. After integrating both curves, we compute the maximum and 
minimum values for the two variables p and A reached during the integration. In Fig. [2l we show 
the positions of those quantities, over the surfaces of section defined in subsection 13.11 With 
those values, first we obtain the positions of the centers for the two orbits through 


Cavrg (C'(A,avrg)) a,vrg)) 


(Max Aavrg Min Aavrg) (Max Pavrg Min Pavrg) 


(19) 


and 


Cnum (C'(A,num)) C'(p^num)) 


(Max Ajium AlinAnum) (Alaxpnum Alinpjmni) 


( 20 ) 








Trojans orbits reproduced by an integrable Hamiltonian 


13 


0.015 

0.01 

0.005 

0 

-0.005 


- 0.01 

-0.015 

- 0.02 

-0.025 


Figure 2: Location of the maximum and miminum values for the coordinates p and A in the 
integration of the fictitious initial condition associated to 1872 Helenos, around L5. Red (light 
gray) points correspond to the averaged level curve and blue (dark gray) points to the numerical 
surface of section 
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The dispersion between centers, which shows how much displacement there is between the orbits, 
is evaluated by putting 


5C 


(C'(p,num) CI(pa,vrg)) ^ C^Ajavrg)) 


c 


(pjnum) 


c 


(A,num) 


( 21 ) 


Furthermore, we obtain the value of the enclosed areas for each curve, as follows: we first 
refer all the points composing the curve, to the already computed center, by introducing the 
quantities 5\ = X — C\ and 5p = p — Cp] therefore, we obtain the distance to the center 
d = \/SfP+S^ and the angle 9 with respect to the horizontal line {9 = atan 2 ((ip, (5A)). Re¬ 
ordering the points by increasing value for the angle 9, we compute the area contained in the 
triangle generated by two consecutive points and the center. The sum along all the triangles 
represents the contained area within each curve, ^num for the complete system and Aavrg for 
average Hamiltonian flow. For a further comparison, we compute also the relative difference 
between the areas JHMn um = l^num — ^avrgl/^numj and the displacement of the centers of 
the orbits with respect to the position of the equilateral Lagrangian point they librate around 
{p = 0, Al 4 = 7 r /3 for L 4 and Xl 5 = Svr/S for L 5 ). 

For 6 of the 7 cases mentioned above, we are able to obtain averaged orbits that reflects the 
behavior of the numerical integrations. In Tabled! where we report the results for the previously 
defined quantities, we show that the averaged areas clearly match their associated numerical 
areas, with a relative error smaller than 2%, except for one case (asteroid 2146 Stentor), for 
which the error is about 13%. This may be due to the fact that 2146 Stentor presents the 
largest displacement with respect to the corresponding Lagrangian point, in a quite anomalous 
orbit. For the rest of the asteroids, the position of the orbits in the surface of section turns 
to be very close to that generated by the equivalent averaged level curve, both centered at a 
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Table 1: Summary of the results for the quantities defining each averaged and numerical or 


ait 


Asteroid 

A 

-^num 

(jA/Anum 

50 

'^(p,num) 

CfA.num) -^1/4,L5 

1868 

2.03 X 10"^ 

3.21 X 10"^ 

6.95 X 10"^ 

-1.08 X 10"^ 

-0.163 

1872 

3.75 X 10-^ 

1.39 X 10"^ 

5.14 X 10"^ 

-6.86 X 10"^ 

-0.235 

2146 

1.67 X 10"^ 

1.25 X 10"^ 

3.71 X 10"^ 

-1.94 X 10"^ 

-0.530 

2207 

2.31 X 10"^ 

6.59 X 10"^ 

7.50 X 10"^ 

-1.31 X 10"^ 

-0.196 

2674 

3.56 X 10"^ 

1.51 X 10"^ 

3.61 X 10"^ 

-1.43 X 10"^ 

-0.077 

2759 

2.67 X 10"^ 

1.29 X 10"^ 

1.04 X 10"^ 

-1.63 X 10"^ 

-0.232 


triangular Lagrangian point. On the other hand, in Table 1 we do not present data for the 
highly inclinec0 (39°) asteroid 2363 Cebriones of our sample. For this asteroid, our normal 
form fails to provide an accurate orbit, using the initial conditions provided in uni. However, 
we find that the numerical orbit generated by 2363 Cebriones presents a very peculiar angular 
excursion (in A) with respect to the Lagrangian point. This failure of the normal form could be 
generated in the initial condition by a non-consistent rotation to the plane of the primaries, in 
the original work. 

4 Summary and perspectives 

In this paper we present a novel normalization scheme, that provides an integrable approximation 
of the dynamics of a Trojan asteroid Hamiltonian in the framework of the Planar Circular 
Restricted Three-Body Problem (PCRTBP). This new algorithm is based on three co-related 
points: the introduction of a set of variables which respects the physical configuration of the 
system; the existence of two degrees of freedom with well distinguishable roles, one corresponding 
to a fast motion and another corresponding to a slow motion, that allows us to fruitfully average 
over the fast angle; finally, the analytic singularity of this model exclusively related with the 
slow angle. 

These three concepts motivates a new way to deal with the initial expansion of the Hamilto¬ 
nian, as it is necessary for the normalization procedure. The slow angle A does not affect the 
solution of the homological equation, determined in order to remove the dependence on the 
fast angle. Thus, we are able to carefully approximate level curves that represent tadpole and 
horseshoe orbits by keeping, in the expansions, a non polynomial dependence just with respect 
to A. 

In order to examinate the accuracy of the normal form produced, we develop some tests. We 
study numerically integrated surfaces of section, along the flow of the complete Hamiltonian, 
and we contrast them with the level curves provided by our integrable normal form, with very 
good agreement. Furthermore, we estimate some quasi-integrals of motion (by computation of 
enclosed areas), which also show excellent agreement with those corresponding to the numerically 
computed surfaces of section. On the whole, this novel approach for producing a new normal 
form results in a very promising approximation of the global behavior of the Trojans motion. 

From a rather theoretical point of view, we think that our normalizing scheme can be com¬ 
plemented with a scheme of estimates, so as to make the remainder exponentially small on a 
suitable open domain. If our algorithm is joint with a Nekhoroshev-like approach, it should 
ensure that the eventual diffusion is effectively bounded (i.e., for intervals of time comparable 

According to the Bowell Catalogue at 2452600.5 JD 
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with the age of our Solar system), for a set of initial conditions. We expect such a set to be 
significantly larger than those considered in the works already existing in the literature. This 
expectation is due to the fact that our method offers a wide coverage of the Trojans orbits. 
For the same reason, we think that our integrable approximation can be used in many cases to 
successfully start the Kolmogorov normalization algorithm. While the corresponding solution is 
valid for any time, the construction of an invariant KAM torus is an extremely local procedure, 
because it must be adapted to the orbit of each Trojan body to be studied. In this context, 
we think that our algorithm can be efficiently used jointly with a KAM-like approach in those 
cases for which the previous implementation of the Kolmogorov normalization algorithm failed. 
Let us remark that such a new application of the KAM theory would require to preliminarly 
construct the action-angle coordinates also for the slow degree of freedom, before starting the 
final normalization procedure. As an alternative strategy, a different formulation of the KAM 
theorem that is not strictly based on action-angle variables could be used [5]. 

More practically, in our opinion the most exciting point is that our method is suitable to 
be translated to more complex models, without requiring essential changes. Since the PCRTBP 
corresponds to a very simplistic representation for the Trojans domain, we are presently working 
to extend the normalization to a Hamiltonian that also considers the eccentricity of the primary. 
The first preliminary results show that our method can be used so as to locate the main secondary 
resonances within the 1:1 MMR region. In particular, we find a good agreement with other purely 
numerical indicators, when the mass ratio between the primaries is small. We plan to include 
these and other results in a future publication. Furthermore, we think that our present and 
future contributions will help to hll that still existing gap between the semi-analytic studies and 
the more complete numerical experiments of the stability region. 
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A Details about the formal algorithm constructing the normal 
form 


As discussed in subsection 12.21 the normalization algorithm dehnes a sequence of Hamiltonians. 
This is done by an iterative procedure; let us describe the basic step which introduces 
starting from when both the values of the indexes ri and r 2 are positive. We assume 

that the expansions of that 


.^2 

T 

i>4 

ri —1 / R 2 

+ '[p 

s=l \l=0 
r 2 — l 


+ E4“’(p — 


2 

2 I ^2 


is)f^ ( +n 

2 


1 >R 2 


s r{ri,r 2 -l;s) I 


(n)/ 


1=0 




l 2 

s r{ri,r2-l-,s) 


l>r2 




( 22 ) 




s>ri l >0 


where ^ E Vi^o V 1 > 4, E Pz,* V 0 < 1 < R2 , 1 < s < ri, E Pz,ri V 0 < / < rs , 
j{ri,r 2 -i,ri) ^ V 1 > r2 , g V / > i?2 , 1 < s < Ti and V 1 > 0, s > ri . Let us 

recall that the Hamiltonian = TL written in ([6]) is suitable for starting the procedure with 

ri = r2 = 1. In formula (I22p . one can distinguish the normal form terms from the perturbing 
part; the latter depends on (^, rf) in a generic way, while in the terms of Z type appearing 
in (122]) the fast variables can be replaced by the action F = (^^ -|- 7/^)/2 . The (ri, r2)-th step of 
the algorithm formally defines the new Hamiltonian 77(^1 >^2) applying the Lie series operator 
exp£^^^^(ri) to the previous Hamiltonian -g prescribed in formula (fTOjl . The 

new generating function iT'^Xr 2 is determined by solving the following homological equation 
with respect to the unknown = Xr^ 2 ^\p^C^ 

r , ^ 4- _ yCi) 

(’’ 1)^2 Jr2 ~ ^r2 ’ 


( 23 ) 
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where we require that is the new term in normal form, i.e. = Zr^^'^ (/o, + r/^)/2, A). 

Proposition A.l If z!^^ = + r/^)/2 and G 'Pr 2 ,ri , then there exists a gener¬ 

ating function ^ 'Pr 2 ,ri ® normal form term Zrl^'^ G Vr 2 ,ri satisfying the homological 
equation (ESD. 


We limit ourselves to just sketch the procedure that can be followed so as to explicitly determine a 
solution of ([23]l and, therefore, prove the statement above. First, we replace the fast coordinates 
(?) v) with the pair of complex conjugate canonical variables {z, iz) such that ? = {z — 'z)ly /2 and 
rj = {z + z)/-\/2. Moreover, the homological equation (f23|) has to be expanded in Taylor series 
with respect to {z,iz), using the slow coordinates (p, A) as fixed parameters (because they are 
not affected by the Poisson bracket L (r,)Zi^\ since zi^^ do not depend on them). Therefore, 

we solve term-by-term the equation m m the unknown coefficient® and 

Cmi,m2,m2,/ci,fc2J SUch. that 

xil^\p,z,\,iz) = Y1 E (cos Aj^Hsm A)^^ {P{X)y 

2mi+m2 A:i+fc2<^+4ri—3 
+ms=l j<2/-|-7ri—6 


and 


Ziy\p,z,X,iz) = Y Cmi,m2,m2M,k2,jp""Hz-iz)^HcosXy^{smXy^ {/ 3 {X)y . 

2miH- A:i+fc2<^+4ri—3 
27n2=l j<2/+7ri-6 


At last, we express the expansions above by replacing {zfiz) with (?, p), and we obtain the final 
solutions in the form = Xr7^(P) ?> "^) P) ziff^ = {Pi (?^ + ^^)/2) A). 

The following property of the Poisson brackets is very useful for our purposes; however, since 
its proof just requires long but basically trivial calculations, it is omitted. 


Proposition A.2 Let f and g be two generic functions such that f G Vr,s o,nd g G Vr'^s' > then 
if r -hr' > 2 ^ {/,£/} G 'Pr+r'-2,s+s' , else ^ {f^g} = 0 . 


In order to provide an algorithm easy to translate in a programming language, we are going to 
give explicit formulas for the new terms of type / appearing in expansion ([9]) of Hamiltonian 
}j(ri,r 2 )_ yg introduce the recursive operation a b, where the previously defined quantity 

a is redefined as a = a b. Therefore, we initially define 

j(ri,r 2 ,s) _ jy^X 2 i,s) V / > i?2 when 1 < s < ri or 

V / > r2 when s = ri or V / > 0, s > ri . 


emphasize that the most celebrated problem concerning the convergence of the normal forms, i.e., the 
occurrence of “small divisors”, does not affect our scheme, because the main integrable term in the homological 
equation (I23j, i.e. Zy\ depends just on the fast action -|- rf‘)l2. 
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Then, we consider the contribution of the terms generated by the Lie series applied to each 
function belonging to the normal form part as follows: 


1 


Ari,r2-,s+jri) ^ 7 

h+j(r2-2) ^ j! ' 


(*) 


V 1 < j < j/ , 0<l <lf, 0 <s <ri , 


(24) 


where the upper limits jj and If on the indexes j and I, respectively, are such that 
jf = I + 1 if r2 = 1 , jf = +00 if r2 > 2 , 

If = +00 if s = 0 , If = R2 + I if 1 < s < ri , If = r2 if s = n , 

= -R 2 + 1 if 1 < s < ri , k = r 2 if s = ri , k = 0 if s > ri . 


(25) 


For what concerns the contributions given by the perturbing terms making part of the expansion 
of ([2^ . we have 


f{ri,r2-,s+jri) i_ cl Ari,r2-l;s) 

Jl+j(r2-2) j\ 


V 1 < j < j/ , l>li, S>1 


(26) 


where the limiting values for the indexes, that are jf and b , are dehned in ()25ll . 

The redefinition rules (I24p and (I26p are set so that the new perturbing part generated by the 
Lie series in (|lUp is coherently split in different terms according to their order of magnitude in p 
and their total polynomial degree in the actions. In fact, by applying repeatedly proposition I A. 21 
to the redefinitions in (??)-(l26]), it is possible to inductively verify that V / > 

li, s > 1. Therefore, the terms making part of the Hamiltonian in the expansion ([9]) 

share the same properties with those appearing in (j22p : this ensures that the normalization 
algorithm can be iterated so as to construct ff{ri,r 2 + 2 )^ ^ simple 

prescription _ jf(ri,R 2 }^ V 1 < ri < , allows to deal with generating functions of the 

next order of magnitude, i.e., As discussed in subsection 12.21 by applying repeatedly 

the rules of the generic normalization step, the algorithm finally ends, by determining the last 
Hamiltonian 




